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TRANSFORMATION OF TWO AND THREE-DIMENSIONAL 


REGIONS BY ELLIPTIC SYSTEMS 

The major effort during this contract period has been the analysis 
of finite-difference methods for composite grids. It was observed that 
linear interpolation between grids would suffice only where low order 
accuracy was required. In the context of fluid flow, this would be in 
regions where the flow was essentially free stream. Higher order inter- 
polation schemes have also been investigated. The well-known quadratic 
and cubic interpolating polynomials would increase the formal accuracy 
of the overall numerical algorithm. However, it can also be shown that 
the stability of the algorithm may be adversely affected. Further numerical 
results are needed in order to assess the nature of this instability 
induced by the interpolation procedure. A complete report on composite 
grid schemes will be presented at the Conference on Large Scale Scientific 
Computation. A copy of that paper is attached to this report. One aspect 
of our work which is not discussed deals with the technical procedures in 
implementing interpolation schemes used on composite grids. Currently 
available software is not designed for repeated interpolation at the same 
points. Therefore, in order to maximize the efficiency of our programs, 
the location of the points used in *he interpolation and the coefficients 
in the interpolation formula are computed only once. By storing these 
values, the cost of computing an interpolated value is comparable to the 
cost of applying the difference equation at an interior grid point. 

The recent appearance of the paper by Hoffman (see attached reports) 
and personal communications with other researchers has motivated us to 
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take a closer look at the order of a numerical algorithm on a curvilinear 
coordinate system. The order of a method is, by definition, dependent 
on the manner in which the grid spacing is decreased. The grid spacing 
may be decreased by adding grid points or by moving existing grid points. 

A detailed analysis of order, including many commonly used mapping 
functions, appears in a paper to be presented at the ASME Applied Mechanics, 
Bioengineering, and Fluids Engineering Conference. A copy of that paper 
is also included in this report. 



To be presented at the Conference on Large Scale Scienti fi c Computation , 
Madison, May 1983 (to be published by Academic PressJ^ 


ERROR ANALYSIS AND DIFFERENCE EQUATIONS ON CURVILINEAR COORDINATE SYSTEMS 

C. Wayne Mastin 


1. INTRODUCTION . 

A computational grid must be constructed when solving partial differ- 
ential equations by finite-difference or finite element methods. Presently 
there are many grid generation algorithms. The choice of algorithm will 
depend on the users desire for control over properties such is orthogonality 
of coordinate lines, location of gr : d points, and smoothness of grid point 
distribution. All of these properties may affect the accuracy of the 
numerical solution. A survey of grid generation techniques may be found 
in the article by Thompson et. al [7]. This report will deal with methods 
for deriving difference equations on curvilinear coordinate systems and 
the effect of coordinate systems on the solution. Recent contributions 
dealing with the effect of the grid on truncation error for one-dimensional 
problems have been made by Hoffman [3] and Vinokur [8]. 

The motivation for this investigation can be seen by considering some 
current problems in computational fluid dynamics. Available computers can 
be used to model the flow about a wing-fuselage configuration. When 
additional components, such as fins, stores, and nacelles are added to the 
aircraft, the computational region becomes increasingly complicated. The 
grid can be extremely distorted and special difference formulas may be 
needed due to irregular neighborhood structures as encountered by by Lee 
et . al . [4] and Roberts [5]. In an attempt to limit grid distortion, 
overlapping grids have also been used for complicated regions. Each 
component is endowed with its own local coordinate system, and interpola- 
tion is used in the solution algorithm. Various interpolation procedures 
have been used by Atta [1], Atta and Vadyak [2], and Starius [6]. The 



possible impact of the interpolation procedure on the solution algorithm 
will be investigated. It is noted that the interpolation technique may 
effect the local truncation error as well as the convergence rate of 
iterative algorithms and the stability of explicit algorithms. 

2. FINITE-DIFFERENCE EQUATIONS . 

A curvilinear coordinate system in the xy-plane is understood to be 
the image of a rectangular Cartesian coordinate system in a £,n-plane. The 
induced grid is therefore composed of quadrilateral cells and difference 
equations may be derived by transforming the partial differential equation 



A typical grid cell is indicated : n Figure 1. This derivation gives no 
information on local truncation error, so an alternate derivation will be 
presented. Regardless of the derivation, the difference equations will 
involve derivatives of x and y with respect to £ and n- Since the grid 
may be given only as a set at data points, it will be assumed that these 
derivatives are approximated using differences. It can be shown that the 
exact computation of these derivatives does not increase the accuracy of 
the method. 

Second order central differences are commonly used in the numerical 
solution of partial differential equations. The truncation error depends 
on the grid spacing in the curvilinear coordinate system, and therefore, 
the corresponding grid spacing in the £n-plane will at present be assumed 
unity. The difference approximations with respect to £; and n are then 
= (fU+l,n) - fU-l,n))/2 
f n = (f(5,n+l) - f(e,n-D)/2 
f u = fU+Vn) + f U-l »n) - 2f(e,n) 

f^ = (fU+l,n+l)+f(5-l.n-l) - fU-l,n+l) - fU+1 ,ti-1 ) )/4 

f n = f(C.n+l) + f(?,n-l) - 2f(t,n). 
nn 
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The local truncation error in using these differences to approximate the 
derivatives with respect to x and y is revealed by examining a series 
expansion of the above differences at (xU,n), yU,n)). First derivative 
approximations are much simpler, and they will be considered first. After 
a little algebra, the difference expression f can be represented as 
if If i a 2 f i 

f r = x r 8x + y^y + 2 + 2*V« + y s x s^ 3x3y 
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The terms in this expansion can be separated into three categories. The 

first order terms are used in deriving the difference equations. The 

second order terms are due to the nonuniform spacing and curvature of the 

curvilinear coordinate systems. The remaining higher order terms (HOT) 

are proportional to the third power of the grid spacing, and terms of this 

order would appear even if a uniform rectangular grid were used. Clearly 

the same remarks can be made about the series expansion for f . From 

n 

these comments it follows that one condition for the derived difference 
approximations to be second order, in the sense that the local truncation 
error is the same order as the square of the grid spacing, is the condi- 
tion that the following quotients 


' r *l 
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(3) 


remain bounded as |r J + |r^| approaches zero. A second condition arises 

when examining the form of the truncation error in the approximation of 
3f 3f 9f 

9x and 9y. From (2) it is seen that the truncation error for 9x can be 
written 


Tx = i(y Te - y T n ) 
Jn 5 
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where 


J = 


x r y 

5 n 


• Vs 


and T§ and Tn are the second and higher order terms in (2) and the 
analogous expansion for f . Certainly some lower bound on the rate at 
which J approaches zero is necessary. Let e be an approximation of the 



angle between the coordinate lines measured by 


y r y n 

e = arctan (p-) - arctan [■—) . 

C n 

If the degree of nonorthogonality is limited by the condition 
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I cot el < M, 


then 


,2 1 

J > 


— M + 1 I'c 
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Once it is noted that each term in TS has a factor of either or y^., it 

follows that when the quotients in (3) and cot e remain bounded as 
| r | + |r | approaches zero, the order of the difference approximations 

S T"| 

. jr first order derivatives is preserved on the curvilinear coordinate 
system. 

The truncation error analysis is considerably more complicated for 
second order derivatives. The major conclusions will be derived without 
going into all the technical details. The series expansions for the second 
order differences are given. 

f = + HOT 


where 
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The mixed derivative approximation involves diagonal neighbors and it is 
convenient to introduce the differences 

f s = (fU+l.n+1) - f({-1.n-l))/2v? 

f t * (f({-l,n+l) - f(E+l.n-l))/2>? 
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f ss ■ (f(5+l,n+l) + f(t-l.n-l) - 2f(c.n))/2 
f tt - (f(f-l .n+1) + f(c+l,n-l) - 2f(t.n))/2. 


The series expansion has the form 

f = T^ + T^ + HOT 
£n 


where 


s 2 f 


a 2 f 


T (1) = x 5n ax + y^y + x^sx 2 ♦ (x^ + x^Jaxsy 

a 2 f 
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Once again the terms have been categorized into those used in the deri- 
vation of the difference equations, and T^, the second and third 

( 2 ) ( 2 ) 

order terms due to the curvilinear coordinates, S ' and T v , and the 

higher order terms which are representative of truncation error that 

would be present on a uniform rectangular grid. The two third order terms 
( 2 ) 

in T' ' which are omitted would be obtained by interchanging x and y in 
the two third order terms that are present. Now the differences with 
respect to s and t can be bounded by differences with respect to £ and n. 
Therefore, the only additional conditions, other than those required for 
first derivatives, in order that the truncation error for second derivat- 
ives be the same order as the square of the grid spacing is that the 
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remain bounded as |rj + |r | approaches zero. 

In the numerical solution of a partial differential equation, the 
truncation error may be decreased at a point by adding grid points or 
moving existing grid points. The technique used to decrease the grid 
spacing has a significant effect on the local truncation error as the 
above analysis indicates. This fact is further illustrated in the follow- 
ing one-dimensional examples. Let xU) be an arbitrary fixed grid point. 
Consider a sequence of grids where the distance between x^U) = x(s) and 
its neighbors, x^U-l) and x^U+l), is decreased by a factor of 2 at 
each step. Then 

x W x ^ x (0) 


and a reduction in the order of the local truncation error would occur 
unless the original grid was uniform; i.e., x = 0. Now consider the 
case where the grid is defined by a mapping xU) = f(c). Suppose the 
neighbors of x(s) are x(£-l) = f(;-h) and xU+1) = f(;+h). Then if 
f' U) f 0 and f"(;) exists 


. f"(c) 

h "° x £ 2 [f'(;)] 2 

is finite. In this case the local truncation error is proportional to the 
square of the grid spacing or 0(h ). It has been assumed that the func- 
tion f and the image of x(s), which is denoted by c, are fixed. This 
conclusion, that the order is preserved, would not necessarily hold if the 
function f changed as h+0. 

It was noted above that the degree of skewness in a nonorthogonal 
coordinate system must be limited to maitain the order of the numerical 
algorithm. The effect of skewness can be further clarified by noting that 



Therefore, for first derivative approximations, the local truncaticn error 

varies inversely as the sine of the angle between the coordinate lines. 

It can also be shown that, for second derivative approximations, an in- 

_2 

crease in truncation error by a factor of sin e is possible. If the 

skewness is accompanied by large variations in grid spacing, this factor 
-3 

increases to sin e. The general conclusion is that a moderate degree of 
skewness has little effect on truncation error. The principal disad- 
vantage in using nonorthogonal coordinates is the added complexity of the 
difference equations. 

Certainly this development does not cover all possible discretizations 
of a partial differential equation. However, similar conclusions hold for 
other commonly used difference approximations. In particular, the same 
conditions for maintaining the order of the difference approximation 
suffice when second order partial derivatives of the form 
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are approximated by 


J { \ [ J (y n f C ' Vn )] C - • Vn>V' 


Extreme distortions in grid cells can have an especially serious 
effect in the finite difference analogs of conservation laws. If the 
partial differential equation 


f x * V s < 6 > 

is approximated by 

l fy n • 9x n>5 + (9 \ • fy £>n = 


then the consistency of the difference equation depends on the difference 
between 



converging to zero where u and v are either 
x and y. 

3. FINITE VOLUME EQUATIONS 
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An alternate method for approximating conservation laws has been 
widely used on curvilinear coordinate systems. Generally referred to as 
the finite volume method, it is derived by integrating the differential 
equation (5) over a grid cell and applying the Gauss Divergence Theorem. 

Let C be an arbitrary grid cell with vertices r( c ,n) . r(s+l,n)> 
rU,n+l)» r(s+l, n+1), where r = (x,y). Let fU + p n + j) denote the 
values of the function f at the centroid of C. In the literature, the 
mean value over the cell is sometimes used rather than the value at the 
centroid. Since the difference in the two values is the same order as the 
truncation error, either definition may be assumed. Integrating equation 
(5) over C gives 


f dy - g dx = s dxdy. 

he h 

A third order quadrature may be derived by using the values of f and g at 
the midpoints of the cell sides. In the usual finite volume formulation, 
this value on the cell side is approximated by the average of the values 
on the two cells having the given side in common. Thus the derived 
difference equation is of the form 

4 

I (yf) . (Ay) . - (ug).(Ax) = As (6) 

i=l 1 1 11 

where A is the area of C and, for example, 

(uf ) 1 = |-(fU -pn + ^r) + fU + |-»n + ^-)) 

( ax) i = xU,n) - xU,n+l). 

The effect of the curvilinear coordinate system can be analyzed by con- 
sidering the difference between yf and the corresponding value of f at the 
midpoint of the common side. The Taylor series expansion of this differ- 
ence is, for example, 

\ (f(C ^.n + ^r) + fU + £tn + -jr)) - f(C,n+£ ) = 

2ax^ ,r ' + '?^ X ££^»n + 2 -) + + y £ »n + ^ + HOT 


( 7 ) 
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x^U.n + y) = x(5+|-.n + ^) + x(C“^,n + ^) - 2x( t *n + -^r) , 

with 

xU.n + y) = ^(x(e,n + l) + xU,n)). 

Thus we again see that the order of the quadrature fotmula is preserved on 

the curvilinear coordinate system if the quotients in (3) remain bounded 

as the grid spacing decreases to zero. Here IrJ and jr I a^e the ais- 

£ n 

tances between the cell centroids. The higher order terms, HOT, in (7) 
are the same order as the square of the grid spacing and hence do not 
decrease the third order accuracy that would exist with the midpoint rule. 

Due to the simplicity of the equation (5) which was considered, 
several aspects of the finite volume method have not been mentioned. In 
practice, s is generally the temporal derivative of some physical quantity. 
When (6) is solved for s, it is apparent that some lower bound must be 
placed on A. This is again accomplished by limiting the nonorthogonality 
and is consistent with results for finite element methods where excessively 
thin elements are to be avoided. It is also noted that, in most finite 
volume methods, the equation (6) is implemented in a two-step algorithm to 
produce second order temporal accuracy. 

Basically then, the finite volume methods also require restrictions on 
the grid to maintain the order of the algorithm. However they can be 
easier to implement in cases where many rectangular grids are patched 
together to fill a complicated region. In such cases it is not uncommon 
for a grid point to have more or less than four neighbors. This causes 
no problem in deriving finite volume equations, but special difference 
formulas must be derived when using the finite difference methods as 
described above. The same comment can be made for cases where triangular 
grids are produced by singularities in the curvilinear coordinate system. 

4. PATCHED COORDINATE SYSTEMS 

Thus far we have considered the curvilinear coordinate system at each 
point to be topologically equivalent to a rectangular cartesian coordinate 
system in the plane. As has been shown, a loss of accuracy in the numeri- 
cal algorithm may occur if the grid is severely distorted. Therefore, 
when the region is too complicated, it may be advisable to partition the 
region and construct a separate curvilinear coordinate system for each 
subregion. Most grid generation techniques are flexible enough to permit 


the smooth continuation of coordinate lines from one subregion into the 
next. However, at points where the boundaries of r :veral subregions 
intersect, there may be greater than or less than four neighboring points. 
An inspection will reveal four grid points which have five neighbors In 
the grid of Figure 2. Special techniques must be applied to derive 
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Figure 2. 

difference equations at these points. One simple technique would be to 
select five nearby points and compute a Taylor series truncated after the 
second order terms. This would give a system of five equations which 
would be solved to obtain difference approximations for the first and 
second order derivatives. Unfortunately, the difference equation derived 
by this method would not resemble the difference equations at other points, 
which it is assumed would be derived by the usual change of variables 
formulas. Thus the effect of tnis differencing technique on numerical 
properties such as stability and iterative convergence would be uncertain. 

A second method of dealing with special points caused by grid patch- 
ing is more compatible with the us:ial differencing techniques. Basically 
it involves selecting nearby points to form a local curvilinear coordinate 
system. Let r(£,n) = (x(s.n). yU.n)) be a grid point with an excess or 
deficiency of neighbors. Then four grid points, denoted by rU±l,n) and 
r(c. n±l). are chosen to define the two coordinate directions through 
r(c.n). Four additional points, denoted by r(c ± 1 , n± 1) and r($ ± 1, n+ 1)» 
are chosen from the four quadrants of the new curvilinear coordinate 
system. The nine points to be used in deriving the difference equation 
at rU,n) have been defined. But the coordinate lines and grid cells may 



be far from that of a uniform rectangular grid. Therefore, as has been 
discussed in Section 2, a loss of accuracy is to be expected when the 
usual finite difference equations are employed. In fact, the local trun- 
cation error for the first derivatives will be the same order as the grid 
spacing. Convergence of second derivative approximations cannot ue 
guaranteed as the grid spacing decreases to zero. Desoite this discourag- 
ing note, accurate results have been computed using this technique. The 
inconsistency of the difference approximation for second order equations 
motivated the search for a higher order approximation. A system of five 
equations in the five partial derivatives of the function f can be con- 
structed by truncating the Taylor series expansions of the central differ- 
ences in (1) after the second order terms. The system is written below 
with the notation of Section 2. 
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For a uniform rectangular grid, the usual difference approximations are 
produced. Although it cannot be guaranteed that the coefficient matrix 
is well-conditioned or at least nonsingular, this is suggested by the fact 
that the system (8) is a perturbation of the nonsingular system of equations 
which produces the usu: 1 difference equations. The later system has a 

coefficient matrix whose determinant is J . This technique generates a ' 
nine-point difference equation using the same differences on the local 
coordinate system that are used at the other points of the grid. 



The local truncation error for first derivative approximations is the 
same order as the square of the grid spacing, whereas, the local trunca- 
tion error for second derivative approximations is the order of the grid 
spacing. This result is valid regardless of the coordinate line spacing 
or curvature. It is only assumed that the coefficient matrix is not ill- 
conditioned. 

While it is possible to generate consistent difference equations at 
special points encountered in grid patching, it should be noted that loss 
of accuracy is possible. The condition that grid lines pass smoothly from 
one region into the next also places a restriction on the number and 
location of grid points in each subregion. Coordinate lines which are 
discontinuous, or have discontinuous slopes, at subregion boundaries can 
be used, but this further complicates the problem of deriving accurate 
difference equations. 

5. OV ERLAPPING COORDINATE SYSTEMS 

When dealing with complicated regions, especially multiply connected 
regions, there may be portions of the boundary where a curvilinear coordi- 
nate system can be easily constructed. For example, consider the region 
between a rectangle and a circle. A polar coordinate system is the 
obvious choice near the circle while a cartesian coordinate system would 
be best near the rectangle. Each of these coordinate systems can be 
extended into the region until they overlap and form a covering of the 
region by grid cells r~ illustrated in Figure 3. In general, there may be 
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Figure 3. 


several overlapping grid systems used to cover a particular region. The 

difference equations must couple the solution values on the various grid 

systems. This transmission of information is most frequently accomplished 

by interpolation at those grid boundary points which lie in the interior 

of the region. Several interpolation procedures will be examined along 

with their impact on finite difference methods. 

Let 6^ ^ be a grid with boundary point r n which is contained in some 
(2) u 

grid cell of the grid G v . First, the general interpolation formula 

f(r Q ) = I a f(r ) (9) 

U j=l 3 3 


will be considered where r. denotes a point in G 

3 


( 2 ) 


When the value of f 
at r^ is replaced in a difference equation by the interpolated value, a 
new difference equation results which may have a different local trunca- 
tion error. Conditions on the coefficients a. which will preserve the 

J 

order of the difference equation can be derived by examining the effect on 

a Taylor series at an interior neighbor r of r^. If the value of f at r^ 

is computed from the values at r. by (9), then 
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This series coincides with the actual Taylor series, computed at r Q , 
through first order terms if 

k k k 

£ a = 1, l a.X = X , l a.y. = y Q , (10) 

j=l J j=l 3 J 0 0=1 3 0 


and through second order terms if, in addition. 



I V x / = {x o )2 ' 

j=l J J 


I a.x.y. 
. 1=1 J J J 
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( 11 ) 


It will be assumed that the conditions, indicated in Section 2, for pre- 
serving the order of first and second order difference equations hold 
for the individual grids and G^. In that case, the local trunca- 
tion error for first derivative approximations would be the same order as 
the grid spacing if (10) holds and the same order as the square of the 
grid spacing if both (10) and (11) hold. The order of the local trunca- 
tion error for second order differences would be the same as the grid 

spacing if both (10) and (11) hold. The order would be the square of the 
grid spacing if an additional condition equating the coefficients of third 

order terms in the series was imposed. The condition which equates the 
p th order coefficients can be written as 


I -[x,f[y,] n ' m = [x 0 f[y n ] 


j=l J 


0- 


( 12 ) 


n + m = p, m = 0, 1 , • • • , p. 

Implicit schemes tend to be difficult to implement on regions which 
use several curvilinear coordinate systems. Therefore, most currently 
used algorithms involve the iterative solution of elliptic equations or 
the explicit solution of parabolic or hyperbolic equations. This naturally 
leads one to question the effect of the interpolation equation on itera- 
tive convergence in the first case and on stability in the later case. No 
detailed analysis will be given here, but a few obvious comments are worth 
noting. Diagonal dominance of the coefficient matrix of the difference 
equations is a sufficient, although not necessary, condition for the con- 
vergence of many iterative methods. A system of diagonally dominant 
difference equations will remain diagonally dominant when the interpolation 
equations (9) are appended if 


I WA £i. 
ri 0 


(13) 


However, when this condition is considered with the first equation in (10), 
which is necessary for consistency, it follows that diagonal dominance 
will be preserved whenever 



The stability properties of (9) can also be observed. Suppose the value 
at r Q is computed at t = (n + l)At by 


k 

f(r n ,(n+ l)At) = l a.f(r . .nAt) . 
U j=l J 3 
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(15) 


The von Neumann stability analysis is based on the behavior of an exponent- 
ial solution of the form 

f(x,y,t) = exp(xt) exp(iyx) exp(ivy). 


Substituting in (15) gives the following 
k 

exp(XAt) = J a. exp(iyAx.) exp(ivAy-), 
j=l 3 3 3 

where Ax. = x. - x~ and Ay. = y. - y~. For real y and v, the exponential 

J J * J J * 

solution will remain bounded as n->-°° provided (13) holds. Therefore, 
whenever (14) holds along with (10), the interpolation equations impose 
no additional stability restriction on the numerical algorithm. 

Several different interpolation schemes will be reviewed in light of 
the above remarks. The first scheme is based on the approximation of a 
linear Taylor polynomial. For each boundary point r n of G^, select a 

(?) (2 r 

nearby point r-j of G v . The neighbors of r.| in (j ' are indexed as in 
Figure 4 so that 

^(r-j ) = (f(r 3 ) - f(r 2 ))/2 
f n (r l } = {f(r 5 } " f ^ r 4^ /2 - 



Figure 4. 


These differences can be used to approximate the partial derivatives with 
respect to x and y. If the approximations are substituted in a Taylor 
series, truncated after the linear terms, then the resulting interpola- 
tion formula can be written 


f(r 0 ) = 

u j=1 J 3 

where 


ORIGINS ‘ rnf 

0F POOR Q uW - n f ' 


“2 " -ct 3 «V X 1> y n ‘ (y 0' y l )x n )/2J 
a 4 = -a 5 = ((y 0 - - (x 0 -x 1 )y c )/2J. 

In this case, it easily observed that equations (10) are valid, but (14) 
does not generally hold. Second degree Taylor polynomials have been used, 
but will not be considered here. It is doubtful that they would give 
better results since the approximation of second derivatives from the 
numerical solution may be very inaccurate. There are other interpolation 
schemes for which both (10) and (14) hold and these will be investigated 
next. 

(21 

Let r^ belong to a grid cell C of G ' with vertices which will be 

denoted by r . , j = 1 , 2, 3, 4, as illustrated in Figure 5. There exists 
3 



Figure 5. 

a unique bilinear mapping of the unit square onto the cell C. The mapp 
ing can be given explicitly by 


r = (1 - s)(l - t) r^ + s(l - t) r 2 + t(l - s) r 4 + st r^ t 



where 0 5 s, t <_ 1. If we set r = r^, then the system of two equations, 
in terms of the x and y coordinates, can be solved to determine the solu- 
tion s = Sq, t = t Q . If f is also assumed to be a bilinear function of s 
and t, then 

f(r 0 ) = 0 - s Q )(l - t Q ) f{r } ) + s 0 (l - t Q ) f(r 2 ) + t Q (l - s Q ) f(r 4 ) 
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It is immediately evident that this interpolation formula satisfies both 
(10) and (14). Although equations (11) would not be satisfied, this 
method can be modified to give higher order interpolants. 3asically it 
involves constructing a bicubic mapping of a square grid in the st-plane 
onto the union of the nine cells consisting of C and the eight cells having 
an edge or vertex in common with C. This mapping can be expressed in terms 
of cubic Lagrange interpolating polynomials. The only additional diffi- 
culty is that' the bicubic equations which determine s Q and t Q would now 
have to be solved numerically whereas the values of Sq and tg in (16) can 
be computed exactly. By construction, the coefficients for the bicubic 
interpolation will satisfy (12) for p = 0, 1, 2, 3. Therefore the inter- 
polation scheme would not increase the local truncation error. However, 
(14) would not be valid. 

Interpolation on triangular regions is very popular in finite element 
analysis. Some of those ideas can be adapted to the present problem. 
Suppose each quadrilateral cell is divided into two triangular cells. Then 
r Q will belong to a triangular cell with vertices r^ , r^ t r 3 . Note that 
this would be the case in Figure 5 if the quadrilateral is partitioned by 
the diagonal from r^ to r^. Now the three equations in (10) determine the 
coefficients for the interpolation formula (9) which coincides with linear 
interpolation on the triangular cell. • As long as rg is an interior or 
boundary point of the cell, the condition (14) will be satisfied. The 
accuracy of the interpolation formula can be increased by increasing the 



Figure 6. 


number of interpolation points. Figure 6 indicates the grid points which 
could be used for quadratic and cubic interpolation. In each case, the 
coefficients in (9) can be calculated from the general equations in (12). 
The quadratic interpolant uses six coefficients obtained by setting 
p = 0, 1, 2 and the ten coefficients for the cubic interpolant are found 
by solving (12) with p = 0, 1, 2, 3. A few cautionary notes are in order. 
The linear interpolation polynomial always exists, but for severely dis- 
torted grids, the system in (12) may be singular and the interpolation 
polynomial may not exist in the quadratic and cubic case. Condition (14) 

is also not satisfied in the quadratic and cubic case. 

Several interpolation schemes have been presented for use on over- 
lapping coordinate systems. This does not include all techniques which 
are presently in use. In particular, we have not considered methods which 
interpolate normal derivatives at the boundary points of each grid. There 
is no reason vhy this analysis cannot be extended to cover that case. We 

would only need consider the formula (9) with some of the r. in and 

( 2 ) 1 
the remaining r. in G x . 

6. CONCLUSIONS . 

The procedure for selecting a curvilinear coordinate system must 
necessarily involve a balance of certain requirements. The rate of change 
in coordinate line spacing and degree of skewness should be limited so that 
the formal accuracy of the difference equation is maintained. On the other 
hand, efficient use of grid points mandates the clustering of points in 
regions where the derivatives of the theoretical solution are large. If 
one must use a highly distorted coordinate system or is faced with the 
prospect of connecting many separate curvilinear coordinate systems in 
different subregions, it is generally possible to derive consistent differ- 
ence approximations. While higher order approximations may exist, their 
use may not be necessary or advisable. The grid spacing is only one 
factor in the local truncation error. The other factor is the theoretical 
solution of the partial differential equation. In a region where all 
derivatives of the solution are negligible, the local truncation error will 
be small regardless of the order. Consequently, when solving fluid flow 
problems, the accuracy of the numerical algorithm is most likely to be 
maintained if irregularities in the grid can be confined to regions of 
free stream flow. There are multitudes of examples where the use of higher 
order methods produce inferior results for one reason or the other. In 
connection with the use of interpolation for overlapping coordinate systems, 
it should be recalled that Lagrange interpolating polynomials may be highly 


oscillatory. 

The analysis of error for nonlinear systems of partial differential 
equations solved numerically on large computational grids can never be 
precise. However the quality of a numerical solution can often be judged 
by examining the grid and the point-to-point variation in the numerical 
solution at the grid points and possibly by recomputing the solution on a 
properly refined grid. 
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f ABSTRACT 

} The order of finite difference representations on 

'% general curvilinear coordinate systems is considered in 
,J some detail. It is shown that the uniform grid order 
| is formally preserved on the nonunifortn, nonorthogonal 
| grid in the sense of the error behavior with an increase 
| in the number of points. However, the coefficients in 
I the series expansion may become quite large for some 
| point distributions. Several specific distributions 
^ are evaluated. 

| INTRODUCTION 

i 

4 Numerically-generated, boundary-conforming curvi- 

| linear coordinate systems have now become common in the 
I numerical solution of partial differential equations, 

I allowing very general codes to be constructed which are 
I applicable to regions with arbitrarily-shaped bound- 
,1 aries. Surveys have been given in (1) and (2) , and a 
| source book with both basic exposition and state-of- 
f the-art developments, ( 3) , has recently become avail- 
, able. 

Difference representations on curvilinear coor- 
i dinate systems are constructed by first transforming 
derivatives with respect to cartesian coordinates into 
expressions involving derivatives with respect to the 
l curvilinear coordinate and derivatives of the cartesian 
coordinates with respect to the curvilinear (the metric 
I coefficients). The derivatives with respect to the 
curvilinear coordinates are then replaced with differ- 
ence expressions on the uniform grid in the transformed 
region. 

,4 Considerable attention is appropriately now being 

, focussed on evaluation of the truncation error of dif- 
Q ference expressions on these curvilinear systems, but 
some misunderstandings have arisen regarding the iden- 
tification of the true order of these expressions. The 
"order" of a difference representation refers to the 
! exponential rate of decrease of the truncation error 
| with the point spacing. On a uniform grid this concerns 
f , simply the behavior of the error with a decrease in the 
point spacing. With a nonuniform point distribution, 
there is some ambiguity in the interpretation of order, 
in that the minimum spacing may be decreased either by 


increasing the number of points in the field or by 
changing the distribution of a fixed number of points. 
Both of these could, of course, be done simultaneously, 
or the points could even be moved randomly, but to be 
meaningful the order of a difference representation 
must relate to the error behavior as the point spacing 
is decreased according to some pattern. This is a moot 
point with uniform spacing, but two senses of order on 
a nonuniform grid emerge: the behavior of the error 

as (1) the number of points in the field is increased 
while maintaining the same relative point discribution 
over the field, or (2) the point distribution over the 
field is changed so as to reduce the minimum spacing 
with a fixed number of points in the field. 

On curvilinear coordinate systems, then, the def- 
inition of order of a difference representation is 
integrally tied to point distribution functions. The 
order is determined by the error behavior as the spac- 
ing varies with the points fixed in a certain distribu- 
tion, either by increasing the number of points or by 
changing a parameter in the distribution, net simply 
by consideration of the points used in the difference 
expression as being unrelated to each other. This 
point is essentially what is noted by Hoffman in (A) ■ 
Actually global order is meaningful only in the first 
sense, since as the minimum spacing is reduced with a 
fixed number of points in the field, the spacing some- 
v;here else must certainly increase. This second sense 
of order on a nonuniform grid then is relevant only 
locally in regions where the spacing does in fact de- 
crease as the point distribution is changed. 

The question of order with nonuniform spacing has 
recently been considered by Vinokur (5) , Hoffman (A) , 
and by Thompson (2) . Other studies of error on curvi- 
linear coordinate systems have been reported in (6-7 ) . 
The present diserssion attempts to clarify this ques- 
tion. 

ORDER ON NONUNIFORM SPACING 

A general one-dimensional point distribution 
function can be written in the form 

*(C> - g(|) (0 < < N) 


Cl) 


la Che following analysis, x will be considered Co vary 
from 0 Co 1. Any ocher range of x can be constructed 
simply by multiplying Che distribution functions given 
here by an appropriate constant. With this form for 
Che distribution function, che effects of increasing Che 
number of points in a discretization of Che field can be 
seen explicitly by defining Che values of £ at Che 
points Co be successive integers from 0 Co N. In this 
form, N + 1 is then Che number of points in che dis- 
cretization, so that che dependence of the error ex- 
pressions on Che number of points in Che field will be 
displayed explicitly by N. This form removes the con- 
fusion Chat can arise in interpretation of analyses 
based on a fixed £ interval (0,1) where variation of 
the number of points is represented by variation of che 
interval AC. The form of the distribution fur.ccion, 
i.e., che relative concentration of points in certain 
areas while the total number of points in the field is 
fixed, is varied by changing parameters in the function. 

The transformation of the first derivative is 


given by 


« _£ 


( 2 ) 


if f is approximated by the second-order central dif- 
ference expression we have, since AC * 1 here. 


f E • « f l+l - f i-l> + T e 


( 3 ) 


where T is the truncation error in this difference 
expression, and itl indicates points adjacent to the 
central point, i.e., indicates increments in C. A 
Taylor series expansion In C yields 

n 1 ( —1) ^ 

T - f, - \ l -r f, , + \ l f, , 

5 C n! (n) „ n n! (n) 


n-0 


n*0 
th 


where f . . represents the n C-derivative of f. The 
n ■ 0 ancTn * 1 terms lead to cancellations, so that 
T^ can be written 

T C " "^,( 20 + 1 )! f (2n+l) (4) 

n* i 

Using (3), the difference expression for f on this 
point distribution is 

f - “ (f ... - f . . ) + T (5) 

x 2x^ i+l l-l x 

where now T = — T. is the truncation error in this 
x x C 

difference representation of f^. From (A) we have then 


f , 


x. “ (2n+l) ! (2n+l) 

C n“ i 


( 6 ) 


point distribution. The first ^-derivative follows 
from (2) ; 


f E ' Vx 


(7) 
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Now f^, and Che higher derivatives, depends on C ex- 
plicitly through the £-der ivatives of the metric deriv- 
atives, e.g., x^, and implicitly through the x-depend- 
ence of f. Thus 


d<f,) 3(f ,) 3(f r ) 

2 . — m ( ) + f — ) X 

dC 1 3C 'x 1 3x ; c 5 


or, in operator form 


_d_ _3_ J_ 

dC " 3C X £ 3x 


( 8 ) 


(9) 


For example, 
f “ 

a 


t 3 . 3 x, 

*3£ + X C 3x f C 


, 3 . 3 

( 3C + X C 3x 


)(x C f x } 


x,.f + x.f 
CC x C xx 


In general, then 
f 


(n) 


dlf 

dc n 


, 3 , 3 . n 

( 3t * 5 m 


(10) 


Note here Chat since f has no explicit C-dependence, 
we have 


* x c £ )f 


x.f 

<, X 


as expected. 

The truncation error in f can then be written, 
using (10) in (6), as 

, i a a 2n+l 

T x--f *. (2ibr <£ ♦ x c f (11) 

c n-1 

Note that the binomial theorem cannot be used to ex- 
pand the pojver of the derivative operator here since 
and x — - do not commute, i.e.. 


3C 


while 


C 3x 


3 . 3 3 “ 

( H )(X E 3x> ' x « S * X E T&i 


, 3 . , 3 . 

(X C 3x )( ?C 


c 3C3x 


Here the metric coefficient, x_, is considered to be 
evaluated analytically, and heftce has no error. (The case 
of numerical evaluation of the metric coefficients 
is considered in a later section.) 

Now the series in (6) cannot be truncated without 
further consideration since the ^-derivatives, f._ . , 

are dependent on the point distribution. Thus if 
the point distribution is changed, either through the 
addition of more points or through a change in the form 
of the disribution function, these derivatives will 
change. Since the terms of the series do not contain a 
power of some quantity less than unity, there is no 
indication that the successive terms become progres- 
sively smaller. 

It is thus not meaningful to give the truncation 
error in terms of (--derivatives of f. Rather, it is 
necessary to transform these £- derivatives to x - 
derivatives, which, of course, are not dependent on the 


Thus all permutations of the operator products of 
degree 2n+l will occur in the expansion of the 2n+l 
power of the derivative operator. For example, with 
2n+l * 3, the following eight operator products will 
occur : 

* ( 3£* (x £ (X £ 3x } ’ (X $ 3x } ’ 

( 3£ )<X C ax^sT*’ (x £ 3x )( 3(f )(x £ lx? * 

^ X £ 3x^7^ * (x £ 3x* ( 3? J 


But since f has no explicit ^-dependence, all of these 
operator products having a on the extreme right will 
make no contribution. Therefore, of the above eight- 


products for l!.** 2n+l power term, only four need to be 


considered : 


(x * ax^ ^ax J * ^ax^V^ax* 

Also since there is no explicit {-dependence in f, the 
following relations apply: 


3. 

I , 
C xxx 


«r* xx 

Here m - 1,2,3, and • 1, ■ 3, and a^ * 3, 


*12 “ *22 


» 1, with all the other a^ being zero. 


K*C &* U 


r t * \ m /w 9 , , df 

( { Sx )3f " X (m+1) dx 
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K ™ rUUIl Mw 'series will each contain 2n+l {-differentiations. 
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where x. , , , Indicates the m+1 { -derivative of x. 
(“+!) t m 


Also 


d C 


and, 




*? w - 


5 dt* 


dx 


Wn 


Other combinations appearing can be inferred from these. 

From these relations it follows that(y- + x 
will expand to the sum of products of {-derivatives of 
x, in each of which product the total number of {- 
dif ferentiations is p. All possible combinations of {- 
derivatives appear in these products: 


1 x 


- 2 

V 

x {{ 

- 3 

* 

4 

x a x a x m 

2 

« 4 

x c* 

5 

x aCv x a * 

- 5 

V 

x au x c \a 


2 2 
,X, , X,, X. 


aaz 

The n-term of the series in (11) then is of the form 
. 2n+l m 2n+l a 

72^171 m- Z l r m^ ih X (i) C-1.2....2n*i) <») 


Now from the form of the distribution function (1), 
it Is clear that 

D. 


‘(i) 


N 


(14) 


where the coefficient D. does depend on the point 
distribution function or (1), but not on the number of 
points, N. Therefore, in (12), 


2n+l a , 
n i» 
i-1 X (i) 


2n+l D. a 

iSi ^ 

N 


2n+l i 
i- n i D i 


im 


2n+l . 

i-Vi 


1m 


2n+l 

i-1 im 


,2n+l 


by (14). The truncation error in the difference ex- 


pression for f then is 


n-1 


2n+l 


£, A_ 


d m f 


(2n+l)!N 2n m “ 1 mn dx® 


where the coefficients, A , given by 

sin 

C 2n+l a 
A »n ' i2l °i 


(15) 


(16) 


depend on the distribution function, but not on the 
number of points. The series (15) is thus a power | 

series in the inverse of the number of points in the f 

field. It therefore is possible to truncate the series^ 
as the number of points in the field, N, increases, 5 
with the result 

3 


1 d f 

T - - E. A . “ 

x ,„2 b*1 ml , m 
6N dx 


(17) 


where the a are non-negative integers on the interval 
(0, 2n+l) such that 


2n+l 

.1 , i a 
i-1 im 


2n+.l (m - 1,2, . . .2n+l) 


(13) 


Also 


1,1 


6 i,2n+l » *i,2n+l " ( 2n+1 ^ 6 1>i 


Neither the exponents, a , nor the numerical coef- 
ficients, C^, depend on Cne point distribution. The 
first and last of the C coefficients are unity: 

C, - C,_., - 1. In (12?, x, , is the i ch {-derivative 


£ X. 2 Xs^an example 


o 

term 


;; a # 


0 


for n - 1 we have the 


1 , 9 A 3,3. 
31 ( 3{ ^3x J f 


where, from (16), 


ml 


C 3 t 

57 i."i D i 


im 


(m - 1,2,3) 


and, as noted above, 

C, - C. - 1 , C, - 3 , a 


13 


3 , a 


12 


‘22 


31 


and all other a, are zero. Thus 
im 


A U “ ■ 


'21 


A 31 “ D l 


The truncation error of the difference expression (5) 
can then be written, using (14), as 


T - - — f - — x f - — x 2 f 
x 6 x^ r x 2 x {{xx 6 x { xxx 


(IS) 


which expands to 


Th« first two terms arise from the nonuniform spacing, 
while the last term is the familiar term occurlng with 
uniform spacing as well. 

From (17) it is clear chat Che difference repre- 
sentation (5) is second-order regardless of the form of 
Che point distribution function in Che sense that the 
truncation error goes to zero as l/N^ as Che number or 
points Increases. This means that the error will be 
quartered when the number of points is doubled in the 
same distribution function Thus all difference repre- 
sentations maintain their order on a nonuniform grid 
with any distribution of points in the formal sense of 
the truncation error decreasing as the number of points 
is increased while maintaining the same relative point 
distribution over the field. 

The critical point here is Chat the same relative 
point distribution, i.e., the same distribution func- 
tion is used as Che number of points in the field is 
increased. If this is the case, Chen the error will be 
decreased by a factor that is a power of the inverse of 
the number of points in the field as their number i.s 
increased. Random additions of points will, however, 
not maintain order. This point has also been noted by 
Hoffman in (A) . In a practical vein this means that a 
solution made with twice the number of points as another 
solution will exhibit one-fourth of the error (for 
second-order representations in the transformed plane) 
when the two solutions use the same point distribution 
function- However if the number of points is doubled 
without maintaining the same relative distribution the 
error reduction will not be as great as one-fourth. 

From the standpoint of formal order in this sense, 
then, there is no need for concern over the form of the 
point distribution. However, formal order in this 
sense relates only to t-he behavior of the truncation 
error as the number of points is increased, and the co- 
efficients A in the series (15) may become large as 
che parameters in the distribution are altered to 
reduce the minimum spacing with a given number of 
points in the field. Thus, although the error will be 
reduced by the same order for all point distributions 
as the number of points is increased, certain distribu- 
tions will have smaller error than others with a given 
number of points in Che field, since the coefficients 
in the series, A , while independent of the number of 
ooints, are dependent on Che distribution function. 

Since the numerical coefficients, , in (16) do 
not depend on the distribution function, ?he quantities 
of concern for the n-cerm of the series (15) are 

A 2n+l a . 

r 1 = r A D i “ <■ - l - 2 2n+I) 

m 1 

, 2n+l la . 

1 „ ,, i» a 

- „ — .n. N x Im 

Nx^ i-1 (i) 

2n+l . . 

.1. a. -1 2n+l N 1- x . a im 

• ‘V “ A <«> 


Now for uniform spacing we have D * 0 for i > 2, and 


Chen by (16) , all A are zero ex 
is given by 


^ Vu.n 7 ”" ich 


*2n+l ,n 


C 2n+1 _ *l,2n+l 
- - — - — D. * * 


Thus the contribution to Che truncation error that 
remains with uniform spacing arises from the m * 2n+l 
term of (15). The ratio of the coefficients A to the 
coefficient A_ . , corresponding to the uniform 

spacing error , Is’Bhen , from (16) and (14). 
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A 2n+l,n 


2n+l 1-1 


C .IT, D, 

m 1-1 1 

C„ .. 2n+l a , . 

2n+l i,2n+l 

i-1 i 

g 

a, 2n+l x,,. lm 

0, l " - c n, (-if) 

i m 1*2 l 


Note that this ratio is the ratio of the coefficient 
of d m f/dx m to that of f in the n-term, i.e., the 

(i; n term, of the series (15) for the truncation 
error. The ratios of the terms arising from the non- 
uniform spacing to that from the spacing itself in the 
n-term of the truncation error expansion (15) as a 
power series in the inverse of number of points are 
then 


A, . ,2n+l 

2n+l ,n d f 


. m 2n+l x,., 

- -4* n (-^-) 

'm ,2n+i , i-2 1 i 

d f x 

. 2n+l K 

dx 


Or der with Fixed Number of Points 

The above considerations have been concerned with 
order in the formal sense of the truncation error being 
reduced by a factor equal to a power of N as the number 
of points in the field is increased, wl'.le maintaining 
the same reLative point distt ibut ion. It has been 
shown that all point distributions main^'in formal 
order in this sense, but that some d isr. r ' but ions miy be 
superior to others with a given total number of points 
in the field. Also, comparisons, may be made on the 
basis of the magnitude of the series coefficients, 
ultimately through the quantities given in (20). All 
this was based on a scries expansion of the error in 
ascending inverse powers of the number of points in the 
field, N. 

An alternate sense of order for point distribu- 
tions is based on expansion of the truncation error in 
a series in ascending powers of the spacing, x . This 
can be developed from che series given above ad (15), 
but with D from (14) substituted in the expression for 
A given By (15): 


Now at least one a^ must be greater than or equal to 
unity for each m, and therefore the exponent of Nx. in 
the above expression is noc negative- Since the a^ do 
not depend on che distribution function, we are lead 
by (19) to compare distribution functions on che basis 
of behavior of the following quantities as Che minimum 
value of x on the field goes to zero with fixed N: 


N 1-1 x . 


(i - 1,2,..., 2n+l) 


C 2n+l , a. C 2n+l 0. . a 
, m _ ,„i . ira m - l.i , im 

A mn * ^ i2l (N x (i) ) “ 0 X x (i) ] 


But, by (13) , 


2n+l ia i-l* a la 2n+l 

no ** - D - D 

1-11 1 1 


*&/*«*** ■ : .* • * 


so that 


121 


where Che N 


i-1 


factor was omitted.) 


, 2n+l x,.. i 

a - cdJ" 

am m i-1 i 

t. 


. 2n+l x, i 

* c„<»V JSx <-^> 


im 


(23) 


Then Che series (15) becomes 


evaluation of distribution functions 

As an example of Che application of the. measures 
of order discussed above, ten distribution functions 
were analysed with specified spacing at £ ■ 0. T,.- 

functions and the coefficients discussed above are 
listed in Tables 1-3, using the following notation: 


(X ,2n 2n+l .m 

r . _ ? (&2 i « — 

x n-1 (2n+l) ! m*l mn ^m 


(24) (i) . i-l x (i) 

S = N ~T~ 
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i 
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where 


2n+l x... a. 

» • c ,n (-if) 

mn m '■! i 
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Recall that the numerical coefficients, C , and the 
exponents, a , do not depend on the distribution func- 
tion. However, in contrast to the series of (15), the 
coefficients, B , may be depondenc on the variation of 
the spacing, x ™ n with a fixed number of points. The 
series here is** therefore not a power series in x , and 
cannot be truncated unless the coefficients, B , are 
bounded as the spacing goes to zero with a fixeS number 
of points. A sufficient condition for this is that the 
quantities involved in the ratio of the coefficients 
to that arising with uniform spacing, i.e., (22), 

-Ip- (i - 1,2,...,-) (26) 


be bounded as x. goes to zero with fixed N. Where this 
is the case, the order of the difference representation 
is maintained with the non-uniform point distribution 
in the sense that the truncation error is reduced by a 
factor equal to a power of the spacing as the spacing is 
decreased with a fixed number of points in the field. 

In the specific distribution functions to be con- 
sidered below, it will be seen that it is possible for 
the quantities of (20) to be larger than those of (26), 
:>ut for most functions the reverse is true. The dif- 
leronce between these two approaches Lo order should be 
kept clear. The first approach coi erns the behavior 
of the truncation error as the number of points in the 
field increases with a fixed relative distribution of 
points. The series here is a power series in the in- 
verse of the number of points in the field, and formal 
o. ier is maintained for all point distributions. The 
coefficients in the series may, however, become large 
for some distribution functions as the minimum spacing 
decreases for any given number of points. Evaluation 
of particular distribution functions in this approach 
la based on the quantities of (20). The other approach 
concerns the behavior of the error as the minimum 
spacing decreases with a fixed number of points in Che 
field. Distribution functions satisfying the condi- 
tims (26) maintain order in this second sense and can 
be compared on the basis of these quant ities . This 
second sense of order is thus more stringent. The con- 
ditions of (26) seem to be unattainable, however. 

Conditions equivalent to .:hose given in (20) for 
comparison of distribution functions were also obtained 
by Vinokur in (5) from consideration of appropriate 
length scales in regions of large gradients. (In Chat 
analysis the transformed variable, £, 1> normalized Co 
the interval (0,1) so the number of points in the field 
does not appear explicitly. The correct interpretation 
of the results of (5) with the present form of distri- 
function is the conditions of (20) and not as given in 


with the subscripts 0 and N indicating evaluation at 
£ “ 0 and N, respectively. The table both the 

values of the coefficients, and , at the points 
of minimum spacing, i.e., £ » 0, and the maximum 
values, together with the location of the maximum. 

The relation of these coefficients to the produce NS 
as the minimum spacing S, approaches zero is also given. 
Plots of the coefficients at £ » 0 and the maximum val- 
ues of the coefficients against Che minimum spacing are 
given in Figs. 1 and 2. The variation of the coeffi- 
cients over the field is shown in Figs. 3 and 4. The 
first of these^shows the entire field for a minimum 
spacing of 10” , while the second gives detail of the 
region near thg minimum spacing (£ « 0) for a minimum 
spacing of 10 . The behavior of the coefficients is 

qualitatively the same for different values of the 
minimum spacing. Finally, Fig. 5 shows the variation 
of £ with x, i.e., the point distribution, ch^ entire 
field being shown for minimum spacings of 10 and 10~ , 
while detail of the region near the mir, imum spacing 
(£ “ 0) is shown for minimum spacings of 10" and 10' . 
Here the ordinate, £, can also be interpreted as the 
fraction of the total number of points that fall be- 
tween x • 0 and the local value of x. 

From Fig. 5b it is clear that, of the functions 
considered here, only the exponential, the hyperbolic 
sine, the hyperbolic tangent, and the error function 
are suitable as point distribution functions with very 
small minimum spacing. The quadratic and sine functions 
do not actually achieve the specified spacing of 10 _& , 
and the rest of the functions concentrate essentially 
all of the points at the left boundary. The error 
function gives the smoothest coverage of the field. 

The hyperbolic tangent is nexr in this regard, while 
the exponential and hyperbolic sine give about the same 
distributions in most of the field. Of the four suit- 
able functions the hyperbolic sine concentrates more 
points near the minimum spacing, i.e., the left bound- 
ary. This function also gives the most nearly uniform 
point distribution in the region of high concentration, 
since the second derivative, and hence L(2) a nd l| 2), 
vanishes at £ ■ 0. This vanishing second derivative 
alao occurs with the tangent and arctangent, but these 
functions concentrate too many of the points near the 
left boundary. 

The plots of the coefficients over the field. Fig. 

? and 4, show L^) for the hyperbolic sine rising 
rapidly from zero to quickly level off just above the 
uniform value for the exponential. The hyperbolic tan- 
gent, by contrast, falls from a value close to that at 
which the hyperbolic sine levels off. The error func- 
tion starts a bit higher than the hyperbolic tangent 
but falls faster. All four of these functions give 
essentially uniform values of the coefficient L^^in 


r« 


ZHCJ Olihi la uOb^'- 


*-• : : L ^:i(.;^s 


the regiqn extending 100 times the minimum spacing from 
the left boundary, cf. Fig. 4a, except for the initial 
rise from zero that occurs for the hyperbolic sine In 
the region extending 10 times the minimum spacing from 
the boundary. The value that occurs for the error 
function is about twice that for the others. (Xitside 
this boundary layer region near the left boundary, the 
hyperbolic tangent and the error function drop off to 
zero, while the exponential and hyperbolic sine remain 
uniform. Thus the hyperbolic sine has the best be- 
havior very near the minimum spacing, while the error 
function, followed closely by the hyperbolic tangent, 
behaves best outside the boundary layer region. The 
error function is, however, a bit higher than the 
others within the boundary layer. Note that the ex- 
ponential, although a suitable distribution function, 
maintains the uniform value near that from which the 
hyperbolic tangent drops off, and therefore the 
exponential is never as good as the hyperbolic 
tangent in regard to the coefficient The 

trends for L^^) ate essentially the same as for 
t (2) f except that now the value for the hyperbolic 
sine is uniform, so that this function has no ad- 
vantage in regard to '. 

For the coefficient Lg , all four functions give 
very nearly the same values within the boundary layer, 
except for the rapid initial rise from zero that occurs 
for the hyperbolic sine and a slightly larger initial 
value occurring for the error function. Outside the 
boundary layer the values for the error function and 
the hyperbolic tangent drop off to zero, the drop 
being a bit faster for the error function, while the 
values for the exponential and hyperbolic sine drop off 
together to a nonzero value. Again the behavior of 
L^' is qualitatively the same. 

It thus appears that the following conclusions can 
be reached on the basis of these coefficients: 

(1) The v v oonential is not <is good as the hyper- 
bolic tangent and therefore should not be used. 

(2) The hyperbolic sine is the best function in 
the lower part of the boundary layer. Other- 
wise this function is not as good as the hyper- 
bolic tangent. 

(3) The error function and the hyperbolic tangent 
are the best functions outside the boundary 
layer. Between these two the hyperbolic tan- 
gent is the better within the boundary layer, 
while the error function is the better outside. 

(4) The logarithm, sine, tangent, arctangent, in- 
verse hyperbolic tangent, quadratic, and also 
the inverse hyperbolic sine (not included in 
Table 1 or the figures) are not suitable. 

... Figs 1 and 2 show that the the variations of both 
and with the minimum spacing are essentially 

same for all four of the suitable functions (ex- 
cept that L^2) an d L( 2 )at 5=0 remain zero for the 
hyperbolic sine). These figures also show that consid- 
eration of the values at 5 = 0 only would be deceptive, 
leading incorrectly to preference for the tangent and 
arctangent, both of which are shown by the other fig- 
ures to be unsuitable. Finally, Fig. 2 shows that the 
four suitable functions do in fact preserve order in 
the sense of variation of truncation error as the num- 
ber of points in the field increases, since has 

only small variation with the minimum spacing. The 
same cannot be said, '..owever, for order in the sense 
of variation of the error as the minimum spacing de- 
creases with a fixed number of points. In fact, Flea. 
2(c) and (d) show that the logarithmic slope of ' 
is near -1 for these functions, and hence the order 
is strictly only first in this sense (since 
the are the coefficients of the x 2 term in the 

error expansion). 


7 


S 



Vinokur (5) considered all of the functions in- 
cluded here, except the exponential, logarithm, and 
quadratic, and also considered the arcsine, which was 
found to be unsuitable. As noted above, the analysis 
of that reference is based on the quantities Lp*. 
Vinokur also shows how to use a basic distribution 
function, with specified slope x at one boundary, to 
construct a distribution function that allows the 
slope to be specified at both boundaries. Forms that 
allow the slope to be specified at an interior point 
are also given. 

Although, as has been shown, all distribution func- 
tions maintain order in the formal sense with nonuni- 
form spacing as the number of points in the field is 
increased. The results obtained for these particular 
distribution functions show that considerable error 
can arise with nonuniform spacing in actual applica- 
tions. Recal that the ratio of the coefficients from 
the nonuniform spacing in the series (15) to the co- 
efficient arising from spacing itself is given by (22), 
which with the definitions of ' 'i'~‘- 
bound for this ratio: 


dx” 


d m f 


Lg i; gives the following 
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A- . , . 2n+l _ — 

2n+l,n d f 


2n+l , . im 

n (p 1 '! 
2n+l, i=2 1L S ’ 


dx 


(27) 


dx 


2n+l 


dx 


2n+I 


The n=l term then yields, for the coefficients involved 
in the leading term of the series (15): 


^31 ^xx» 


= lP> 


A,, f 
23 xx 


*31 


= 3L 


( 2 ) 


XXX 


(28) 


Now a typical case involving a boundary layer might 
have 100 points with a minimum spacing of 10 rela- „ 
tive to.a maximum field extent of unity. Thus N = 10^, 
S a 10 ", and NS = 10 . Then for 


(i; 


^NS' 


as for the best of the functions considered, we have 


( 2 ). 


10 


l< 3) . 10 s 


and then the ratios of the error from the nonuniform 
spacing to that which arises from the spacing itself 
are, approximately, 

m 8 f x , . „4 f x 
10 -= and 10 


XXX 


XXX 


n 

Sinjfj the error term from the spacing here is S f 

^ xxx » t *' e err °r terms due to the nonuniform XXX 
spacing are 

■ io- 8 f 


ltT’f 


and 


as compared with 10 _12 f due to the spacing. Now 
for the same number of points with uniform spacing we 
would have a spacing of lO - ^ anc j an error of 10“^f 
Thus the error due to the nonuniform spacing in this* 
case is well below what would occur on a uniform grid 
with the same number of points, except for the f term. 
(It will be shown below that this term can be elimi- 
nated from the truncation error by evaluating the co- 
ordinate derivatives numerically rather than analyti- 
cally.) 

This example shows that the contributions to the 
error from the nonuniform spacing are significant and 
must be considered. While the contribution form the 
spacing Itself decreases with the spacing, the con- 
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Cributions from Che nonuniform spacing increase as the 
spacing decreases for very small spacings. 


The lead term of the error then is 
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All of the above considerations have assumed that 
the derivatives of x with respect to £ are evaluated 
exactly. If the coordinate derivative, x , in the 
difference expression (S) is evaluated numerically by 
the same central difference expression used for f we 
have, in place of (5): 

f - — + T' (29) 

x x i+1 - x ±1 x 

With x expanded in Taylor series we have 

o 4 . ? ? X (2n+1) 

x i+l " x i-l 2x £ n=l (2n+l) ! (30) 

Using this and (5) in (29) we then have 

f - T 

f - - + T* (31) 

X i + ± ? X (2n+1) X 

x^ n“l (2n+l) ! 

Now the f^term of T^ corresponds to m = 1 in (15) : 

I A ln 

“7 1 (2n+l)!N 2n 


ard by CL® , 


C. 2n+l a., 

A. - IT 
In i=l l 


But C, = 1 and a., = 6, _ as given above, so that, 
1,,,, il i,2n+l’ 

usrng (14), 

°2n+l N X (2n+1) 

A in “ 

Then the coefficient of f in T is 

x x 


which is the same as (17), except that the lower limit 
of m is 2 in the present case. This can finally be 
written as 

T' - - 4 x__f - z x 2 f (33) 

x 2 ££ xx 6 £ xxx 

Thus the use of numerical evaluation of the co- 
ordinate derivative, rather than exact analytical 
evaluation, eliminates the f term from the truncation 
error. Since this term is tfie most troublesome part 
of the error, being dependent on the derivative being 
represented, it is clear that numerical evaluation of 
the metric coefficients oy the same difference repre- 
sentation used for the function whose derivative is 
being represented is preferable to exact analytical 
evaluation. It should be understood that there is no 
incentive, per se, for accuracy in the metric coef- 
ficents, since the object is simply to represent a 
discrete solution accurately, not to represent the 
solution on some particular coordinate system. The 
only reason for using any function at all to define the 
point distribution is to ensure a smooth distribution. 
There is no reason that the representat ions of the co- 
ordinate derivatives have to be accurate representa- 
tives of the analytical derivatives of that particular 
distribution function. 

Two-Dimensions 

The two-dimensional transformation of the first 
derivative is given by 

f x ’ X f 5 - y £ V <34) 

where the Jacobian of the transformation is 


J = x.y - x y. 

a n n >, 


_L ? ( 2n+l) 

' x^ n=l (2n+l) ! 


and, for use in (31), we have 


With two-point central difference representations for 
all derivatives, we have 

6y6f - j yj f 

f = J31— 5 i— -0_ + T (35 ) 

x i ? x« n y -S n x 6 c y x 


, 1 ® (2n+lK 

f x ~ T x ■ f x (1 + x^ n=l (2n+l)! ) 

, 2n+l ,m, 

+ f 1 z A ±1 

n " 1 (2n+l)' !N 2n m=2 2111 dx” 


But the coefficient of f on the right here is exactly 
the denominator in ( ) , so that, using (14) in this 

denominator, we have the following expression for the 
truncation error in the difference representation (29): 


(2n4-l) !N 2n m-2 A mn dx m 


1+i l 



D 1 n * 1 (2n+l)!N 2n 


where 


6 £ f 5 f i+l,j " f i-l,j 


V ' f i,j+l ' f i,j-l 


and is the truncation error. After expansion of all 
quantities in Taylor series about the central point 
and considerable algebraic manipulation, we have for 
the leading term of the truncation error 

T x - S<Wnn - Vn x « )f x* 

*S ( >sV ( > nr, - y ES )f yy 

+ £Vn<\n - x £t> * x n y E y rtn ' Vn y EE If xy 

+ second-order terms in the spacing (36) 

where the coordinate derivatives are understood to 


tart. ! 


aubs. -;aer.t oat;-: s 


represent central difference expressions, e.g.. 


kx 

2' X i+l,j 


- x i _l,j ) * x n " 2 (x i,j+l 




T T 

T X " sTn() n -^) ( 3 i ^n COS ^ ^ • sin V° S *n r> (38) 


^Therefore the truncation error, in general, varies in- 
2x ij + x i-i, ORIGINAL PA3E is} versely with the sine of the angle between the coord i- 
oftOQ QUALITY "ate lines. Note that there is also a dependence on 
2x + x ™ the direction of the coordinate lines. To further 

U clarify the effect of nonorthogonality, the following 

example is included. For simplicity, only the trunca- 
the truncation error arise tion error terms arising from nonuniform spacing are 

;. The familiar terms pro- considered, 

spacing occur in addition The contribution from nonorthogonality can be 

isolated by considering the case of skewed parallel 


These contributions to the truncation error arise 
from the nonuniform spacing. The familiar terms pro- 
portional to a power of the spacing occur in addition 
to these terms as noted. 

Sufficient condit ions can be stated for maintaining 
the order of the difference representations. First of 
all, as in the one-dimensional case, the ratios 


lines with x 
below: 


0 as diagramed 
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Here (36) reduces to 


must be bounded as x. 


approach zero. 


uiuoL uc a o a _ , a , Jr* J • » 

second condition must be imposed which limits the rate 
at which the Jacobian approaches zero. This condition 
can be met by simply requiring that cot9 remain bounded, 
where 9 is the angle between the £ and n coordinate 
lines. The fact that this bound on the nonorthogon- 
ality imposes the correct lower bound on the Jacobian 
follows from the fact that 


T x 2 X £S f xx + 


y f 

nn yy 


if 

Since cot0 = — — 


- 

2Sc c ;x 


55 f xy 


this may be written 


T = - 
x 


X_ _f 
« XX 


+ -|(y 


f 

nn yy 


- x..f )cot0 
CS xy 


implies 

,2 1 , 2 2 22 22 2 2 , 

J :ST lx c\ + Vn + V£ + , « y n’- 

With these conditions on the ratios of second to first 
derivatives, and the limit on the nonorthogonality 
satisfied, the order of the first derivative approxi- 
mations is maintained in the sense that the contribu- 
tions to the truncation error arising trom the non- 
uniform spacing will be second-order terms in the grid 
spacing. 

The truncation error terms for second derivatives 
that are introduced when using a curvilinear coordi- 
nate system are very lengthy and involve both second 
and third derivatives of the function f. However, it 
can be shown the same sufficient conditions, together 
with the condition that 


remain bounded, will insure that the order of the 
difference representations is maintained. 

It was noted above that a limit on the nonortho- 
gonality, imposed by (37), is required for maintain- 
ing the order of difference representations. The 
degree to which nonorthogonality effects truncation 
error can be stated more precisely. The truncation 
error for a first derivative f can be written 

T , ■ JVs - W 

where T and T are the truncation errors for the 
difference expPessions f and f . Now all coordinate 
derivatives can be expressed using direction cosines 
of the angles of inclination, <p and $ of the £ and 
q coordinate lines. After some simplification, the 
truncation error has the form 


The first term occurs even on an orthogonal system and 
corresponds to the first term in (33). The last two 
terms arise from the departure from orthogonality. 

For 0 < 45®, these terms are no greater than those 
from the nonuniform spacing. Reasonable departure 
from orthogonality is therefore of little concern when 
the rate-of-change of grid spacing is reasonable. 

Large departure from orthogonality may be more of a 
problem at boundaries, where one-sided difference ex- 
pressions are needed. Therefore, grids should probably 
be made as nearly orthogonal at the boundaries as is 
practical. Note that the contribution from nonortho- 
gonality vanishes on-a skewed uniform grid. 
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Tabic 1 


Coefficients ct Minimum Spacing 



FUNCTION 

NS 


<L (3) ) 

' N '0 

(L (2) > 

' s ; o 

(L (3) ) 

1 S l 0 

Exponential: 

a* 5 - 1 

a 


2 


. a , . 2 

tC - 1) 

' - 1 

a* - 1 

a 

a 


Hyperbolic 

Tangent 

j tanh a(l - £) 

tanh a 

2a 

ainh2a 

2atanh a 

2o 2 (Hanh 2 a * l) 

2slnh 2 o 

f(3tanh 2 a - l)ainh 2 2a 

Hyperbolic 

slnh a\ 

a 


2 


. .2 
■ lnh a 

Sine 

slnh a 

slnh a 

0 


0 

Error 

Function 

, erf a(l - £) 

erf a 

2 

2_ ae~° 
^erf a 

2a 2 

2a 2 (2a 2 - 1) 

2 

ae° trfu 

j(2a 2 - l)(e a erfa) 2 

Tangent 

t4n (0 < a < — ) 

a 


, 2 


, 2 

tan a 2' 

tan a 

0 

2 i 

0 

2tan a 

Arctangent : 

, tan l a(l - £) 

a 

2a 2 

2 2 

2a (3a - 1) 

2ntan 3 a 

2 (3a 2 - 1) (tan* 1 a ) 2 

. -1 
tan a 

(1 ♦ a^)tan 

1 ♦ a 2 

(i ♦ « 2 ) 2 

Sine : 

, *10 a(l - O * 

a 

a tan a 

2 

„ 2 

2 

1 Tin a (0 ! “ i 2 } 

tan a 

Ct 

tan a 

-tan a 

Log 

, fn[l ♦ a(l - Ol 

a 

a 

2o 2 

fn(l + a) 

2[fn(l + a) I 2 

fn(l + a) 

<1 + a)ln(l + a) 

i + a 

a + o) 2 

Inverse 
Hyperbolic : 
Tangent 

tanh /#x , . 

Zi CO c 0 « l) 

tanh' a 

a 


2a 2 


2(tanh' 2 a) 2 

tanh~^a 


0 

Quadratic : 

a£ + (1 - a)l 2 (0 < a < 1) 


2(1 - a) 


2(1 - a) 


a 

a 


2 



a 


Table 2 


<\ax 


Exponential: same (uniform) 


Hyperbolic 

Tangent 


same ( l 


0 ) 


Hyperbolic 

Sine 

: a tanh a (f t * 1) 

Error 
Func t Ion 

i 

• 

/*»» 

• 

o 

Tangent 

: 2a tana (£ - 1) 

Arctangent 

: o (€ - 1 + 

Sine 

: same (£ ■ 0) 

Log 

j a (1 - 1) 

Hyperbol lc 
Tangent 

: (1 - D 

1 - a* 


Quadratic : saae (( • 0) 


Maximum Coefficients 


u4 3) , 

N max 

<L< 2 >) 

S max 

(l< 3 >) 

S max 

same (uniform) 

same (Q * 0) 

some (C * 0) 

same (£ • 0) 

same ( ( ■ 0) 

same (t * 0) 

same (uniform) 

'i slnh a (slnh af. * 1) 

same (£ * 0) 

same (( » 0) 

same (£ • 0) 

same (£ ■ 0) 

2a 2 (3tan 2 a + 1) 

tan a U - 

9 2 - 1 

7- tan a (tan at * ) 

4 /r 

9 2 ,7 . 1 x. 

same 

same 

same (uniform) 

same (£ • 0) 

same (t - 0) 

2a 2 (C - 1) 

same (uniform) 

Same (uniform) 

2a 2 (3a 2 >1) - 

-1 

2at*nh a (C • D 

2(3a 2 + l)(tanh" l a) 2 (C 

(1 - o ) 

same (? - 0) 

same (£ * 0) 

same (f • 0) 


l) 


KOTE: 'Same' Indicates maximum value is same aa value at T * 0* 



Table 1 


Coefficient! a» Minimum Spacing, Approaches Zero 



«P>. 


3M 


BSM 


(L< 2) ) 

S nax 

(L< J> ) 

S ma 

Exponential: 

t "i5 


M 


B 


same 

same 

Hyperbolic 

Tangent 

^Ts 

t s > 2 

i 

NS 

< i -) 2 

'NS 7 

same 

sane 

same 

same 

Hyperbolic 

Sloe 

0 

& 

0 

(X) 2 

W 

* A 

same 

On^) 2 

same 

Error 

Function 


£) 2 

/r i 

2 NS 

4 l NS> 

same 

same 

same 

same 



2 


2 . „ 

2 



2 

Tangenc 

0 

* 

0 

T ( nI> 

e 1 

(i) 2 

8 NS 

» i 

i=i ( ±) 

16 l NS' 


2 


2 NS 

2 NS 

Arctangent : 

2 

1 

2 

NS 

<s )2 

/3 1 
« NS 

_2_ (X) 2 
, 2 W 
2» 

same 

same 


2 , 

2 

2 

2 . . 





Sine : 

V 1 

4 NS 

It 

4 

4 W 

- T 

same 

same 

same 

same 






i 

i , 



log : 

i 

2 

NS 

2 ^> 2 

NS 

2(. NS > 

same 

same 

Inveree 







2 

„ 1,2 

Hyperbolic : 
Tangenc 

0 

2(1 - NS) 

0 

2(tanh* 1 o) 2 

2 

NS 

8(X) 2 

NS 

NS 

*'NS J 

Quadratic: 

z_ 

NS 

0 

2(^) 2 

NS 

0 

same 

sane 

same 

same 

NOTE: 'aaae 

' indicates 

maximum value 1 

s same as value at £ ■ 0 

i. 
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Fig. 5 Point Distribution 
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